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Abstract. By means of three-dimensional hydrodynamic simulations with an Eulerian PPM code we investigate 
the time-dependent evolution and properties of accretion tori around nonrotating and rotating stellar-mass black 
holes, using a pseudo-Newtonian (Paczyriski & Wiita or Artemova-Bjornsson-Novikov) potential to approximate 
the effects of general relativity. The simulations are performed with three nested Cartesian grids to ensure sufficient 
resolution near the central black hole on the one hand and a large computational volume on the other. The black 
hole and torus are considered as the remnant of a binary neutron star or neutron-star black-hole merger. Referring 
to results from previous hydrodynamical simulations of such events, we assume the initial configurations to consist 
of a black hole with a mass of about 4Mq girded by a toroidal accretion disk with a mass in the range from 
about 0.01 M© to 0.2 M®. We simulate the torus evolution without and with physical shear viscosity, employing 
a simple a-model for the gas viscosity. As in our previous work on merging neutron star binaries and neutron 
star/black hole binaries, we use the equation of state of Lattimer and Swesty. The energy loss and lepton number 
change due to neutrino emission from the hot torus are treated by a neutrino-trapping scheme. The energy 
deposition by neutrino-antineutrino annihilation around the disk is evaluated in a post-processing step. The time- 
dependent efficiency of converting gravitational energy to neutrinos, expressed by the ratio of neutrino luminosity 
to accretion rate of rest-mass energy, can reach maximum values of up to about 10%. The corresponding efficiency 
of converting neutrino energy into a pair-photon fireball by neutrino annihilation peaks at values of several percent. 
Interestingly, we find that the rate of neutrino-antineutrino annihilation decays with time much less steeply than 
the total neutrino luminosity does with the decreasing gas mass of the torus, because the ongoing protonization of 
the initially neutron-rich disk matter leads to a rather stable product of neutrino and antineutrino luminosities. 
The neutrino luminosity and total energy release of the torus increase steeply with higher viscosity, larger torus 
mass, and larger black hole spin in corotation with the disk, in particular when the spin parameter is a ^ 0.8. 
The latter dependence is moderated in case of a high disk viscosity. For rotation rates as expected for post-merger 
black holes (a 0.5) and reasonable values of the alpha viscosity of the torus (a ~ 0.1), torus masses in the 
investigated range can release sufficient energy in neutrinos to account for the energetics of the well- localized short 
gamma-ray bursts recently detected by Hete and Swift, if collimation of the ultrarelativistic outflows into about 
1% of the sky is invoked, as predicted by recent hydrodynamic jet simulations. 
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1. Introduction 

Merging double neutron stars or neutron-star black-hole 
binaries have long been discussed as possible sources 
of cosmic gamma-ray bursts (GRBs) (e.g., Blinnikov et 
al. 11351 Paczynski |P9"8l>l Goodman IT^gfjl Goodman et 
al. 115571 Eichler et al. iTP8"9l Paczynski ITOOTl Narayan 
et al. IT5321 Meszaros & Rees 113331 Woosley UQ9^al 
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Jaroszyhski ITMfl [15551 Mochkovitch et al. IT5531 IT5551 
Thompson HHMI Witt et al. lT551 Janka & Ruffert IT5551 
Meszaros & Rees 119971 Popham et al. 119991 Meszaros et 
al. I1999[) . After the discovery that long GRBs are linked 
to massive star explosions, which confirmed predictions by 
Woosley l |1993af and MacFadyen & Woosley l(T555}l . it is 
in particular the class of short, hard GRBs with durations 
of less than about 2 seconds, which is thought to origi- 
nate from colliding compact objects in binary systems. 
Recent observations of well-localized short GRBs (Gehrels 
et al. 1217051 Hjorth et al. 1217051 Berger et al. l2U05jl seem to 
support this hypothesis, because in agreement with expec- 
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tations short GRBs are found in galaxies without strong 
star-formation activity. 

In the most popular scenario, GRBs have in common 
that they signal the birth of stellar-mass black holes and 
the associated huge energy release when the growing black 
hole accretes matter at enormous rates between fractions 
of a solar mass per second up to many solar masses per 
second fWooslev ll993al Popham et al. 119991 MacFadyen 
& Wooslcy 1999). The huge amount of gravitational bind- 
ing energy released during the accretion process of up to 
several solar masses of gas into the black hole could ex- 
plain the energetics of even the most distant cosmologi- 
cal gamma-ray bursts (e.g., GRB981214, see Kulkarni et 
al. I1998[) . Moreover, the compactness of the stellar-mass 
black hole could naturally produce the rapid variability 
on time scales of milliseconds observed in many bursts. 
For these reasons, massive accretion disks or thick accre- 
tion tori around stellar-mass black holes are considered as 
favourable model for the central engines of the cosmolog- 
ical GRBs. 

In this scenario the energy of the ultra-relativistically 
expanding fireball or jet, can be provided by the anni- 
hilation of neutrino-antineutrino pairs (Paczyfiski 119911 
Meszaros & Recs 1997 Wooslev ll993bl Jaroszyflski[1993 
HMol Mochkovitch et al. ITM31 IT39l)|> or possibly by mag- 
netohydrodynamical processes (Blandford & Znaiek H977l 
Meszaros and Rees HM?! Popham et al. H^9l Li l2"0"UUl Lee 
et al. 120001 see also Livio et al. 119991 Rosswog & Ramirez- 
Ruiz 2002; Rosswog et al. 200SJ. In the former case, the 
gravitational binding energy of accreted disk matter is 
tapped, in the latter case (also) the rotational energy of 
the central black hole (BH) could be converted into kinetic 
energy of the outflow. 

In a series of previous papers (Ruffert et al. H99Hl fl997 
Ruffert & Janka 1999 , 2001) it has been shown that the 
neutrino emission associated with the dynamical phase of 
the merging or collision of two neutron stars is not suffi- 
ciently powerful and too short to provide the energy for 
gamma-ray bursts by neutrino-antineutrino annihilation. 
After a transient phase of gravitational stability, the com- 
pact massive remnant of the neutron star merger may 
collapse to a black hole and some matter may remain 
in a hot toroidal accretion disk around the black hole, 
which continues to radiate neutrinos with high luminosi- 
ties (for recent discussions and model calculations, see, 
e.g., Popham et al. lTMfl Rosswog et al. 113591 12UUU1 Ruffert 
& Janka ITMfl Oechslin et al. l2UUl Morrison et al. lSUMl 
Oechslin & Janka l3UUrJ|l . 

Newtonian hydrodynamic simulations including a 
physical equation of state (Ruffert et al.©96. Ruffert & 
Janka 1999) indicate that about 0.1 Mq of matter might 
obtain enough angular momentum during the merging of 
the neutron stars to avoid immediate accretion by the 
black hole. These numbers have recently been confirmed 
by simulations including a treatment of general relativity 
(Oechslin & Janka 2006). The presence of large angular 
momentum allows the matter to form and sustain toroidal 
disk around the black hole for a period of time depend- 



Table 1. Parameters of the torus evolution models. The 
3D simulations were performed with three levels of nested 
Cartesian grids with the number of grid cells per level 
and spatial direction given in column "zone" . Ai ca i is the 
time interval over which the model was evolved. All mod- 
els started with a black hole mass of Mg H = 4.O17M . 
The initial torus masses are given by M^, < a < 1 is 
the initial value of the dimensionless BH spin parameter, 
"pro" means that torus and BH were chosen to rotate 
in the same direction, "ret" means opposite directions of 
rotation, and a denotes the dimensionless disk viscosity 
parameter. 



Model 


zone 




Mi 


BH spin 


direction 


vise 






HIS 


Mr-, 


a 




Ot 


r00-32 


32 


70.0 


0.0478 


0.0 




0.0 


all-32 


32 


70.0 


0.0478 


0.0 




0.001 


al2-32 


32 


70.0 


0.0478 


0.0 




0.004 


al3-32 


32 


70.0 


0.0478 


0.0 


- 


0.01 


al4-32 


32 


70.0 


0.0478 


0.0 


~ 


0.1 


rol-32 


32 


70.0 


0.0478 


0.3 


pro 


0.0 


ro2-32 


32 


70.0 


0.0478 


0.6 


pro 


0.0 


ro3-32 


32 


70.0 


0.0478 


0.8 


pro 


0.0 


ro4-32 


32 


70.0 


0.0478 


0.3 


ret 


0.0 


ro5-32 


32 


70.0 


0.0478 


0.6 


ret 


0.0 


irl-32 


32 


70.0 


0.0120 


0.0 




0.0 


ir2-32 


32 


70.0 


0.0239 


0.0 




0.0 


ir3-32 


32 


70.0 


0.0956 


0.0 




0.0 


ir4-32 


32 


70.0 


0.1912 


0.0 




0.0 


r00-64 


64 


40.0 


0.0478 


0.0 




0.0 


al3-64 


64 


40.0 


0.0478 


0.0 




0.01 


al4-64 


64 


40.0 


0.0478 


0.0 




0.1 


ri4-64 


64 


40.0 


0.1912 


0.0 


pro 


0.0 


ro2-64 


64 


40.0 


0.0478 


0.6 


pro 


0.0 


ro5-64 


64 


40.0 


0.0478 


0.6 


ret 


0.0 


irl-64 


64 


40.0 


0.0120 


0.0 




0.0 


ir4-64 


64 


40.0 


0.1912 


0.0 




0.0 


ir5-64 


64 


40.0 


0.1912 


0.0 




0.1 


arl-64 


64 


40.0 


0.0478 


0.6 


pro 


0.1 


ar2-64 


64 


40.0 


0.1912 


0.6 


pro 


0.1 



ing on the viscous transport of angular momentum in the 
disk. Therefore there could be enough time for this ma- 
terial to radiate away in neutrinos a fair fraction of its 
gravitational binding energy. A similar situation could re- 
sult from the merging of a neutron star with a black hole 
(Lee & KluzniakUigEHjll Kluzniak & Lee lTMSl Janka 
et al-EMi Rosswog et al. l2lMjl . 

Many questions remain to be explored by quantitative 
calculations of the formation, evolution, and energy re- 
lease of accreting stellar-mass black holes as possible cen- 
tral engines of GRBs. The results of the three-dimensional 
(3D) hydrodynamic simulations presented in this paper 
aim at studying the time-dependent evolution of such 
black holes girded by accretion tori that were formed from 
the decompressed matter of a neutron star, which was 
disrupted in a compact binary merger. Such tori are gen- 
erated by a violent dynamical process in the immediate 
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vicinity of the black hole. The assumption of stationary 
accretion is therefore valid only approximately and only 
after a longer phase of adjustment and relaxation of the 
matter girding the black hole. This relaxation can pro- 
ceed over a fair fraction of the whole period of massive ac- 
cretion and requires time-dependent modelling. Therefore 
our simulations intend to contribute to answering the fol- 
lowing questions: How much mass remains in the accre- 
tion torus for longer than dynamical time scales? What 
is the mass accretion rate of the black hole as a function 
of time? What are the properties of the accretion torus, 
its density, temperature, neutron-to-proton ratio, neutrino 
luminosity as functions of time? How do they depend on 
torus mass, gas viscosity, and black hole rotation? How ef- 
ficient is neutrino-antineutrino annihilation in depositing 
energy in the surrounding of the black-hole torus system 
in dependence of these parameters? What are the implica- 
tions for powering ultrarelativistic outflow and gamma-ray 
bursts by neutrino-antineutrino annihilation? 

Our 3D models of the time-dependent, non-stationary, 
evolution of the black- hole torus systems include a detailed 
treatment of the equation of state and neutrino physics. 
They thus replace approximations of previous work as 
made in the analytic studies of Popham et al. (|1999p and 
Di Matteo et al. (20@ and in the numerical 2D mod- 
els of Lee & Ramirez- Ruiz (2002 ). These latter simula- 
tions with azimuthal symmetry employed a simple ideal 
gas equation state and made the assumption that all the 
dissipated energy is radiated away in neutrinos. Lee et 
al. H20041 12005al I2005b|) . used a more realistic equation 
of state and took into account neutrino energy loss and 
opacity effects, but with a simplified model for changes 
lepton number (Y e ), considered only Newtonian gravity 
and nonrotating black holes, and did not evaluate their 
2D hydrodynamic models for the energy deposition by 
neutrino-antineutrino annihilation. 

The paper is organized as follows. In Sect. [5] we will 
introduce the numerical tools used for our simulations, 
outline the initial conditions, and give an overview over 
the different investigated cases. Sect. will present the re- 
sults of our simulations and discuss their implications. The 
summary of our main findings will be given in Sect.0 and 
concluding remarks and an outlook will follow in Sect. El 

2. Numerical Models 

2.1. Dynamics 

The three-dimensional computations of black-hole 
accretion-disk dynamics in this paper were performed 
with a Newtonian hydrodynamics code based on the 
Piecewise Parabolic Method (PPM) of Colella & 



Woodward (1984} with three levels of nested grids 
(Ruffert 1993J) to ensure both sufficient resolution near 
the central black hole and a large computational volume. 
A higher grid level has twice the zone size of the level 
below. With an equatorial length and width of the total 
computational volume of 500 km, the smallest zones 
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Fig. 1. The black hole event horizon th, the innermost 
stable circular orbit nsco , and P as defined in Eq. (j2J , as 
functions of the black hole spin parameter a. Both th and 
nsco are given in units of GMbh/c 2 = 0.5i? s . 

have a side length of 1.95 km for 64 zones per grid and 
dimension ("64- resolution") and 3.90km for 32 zones on 
each grid level ("32-resolution"). 

The "black hole" with mass Mbh is treated as a grav- 
itational centre surrounded by a vacuum sphere on the 
computational grid. Its gravitational potential ^bh as a 
function of radius r is calculated from 



d$ 



BH 



GM] 



BH 



dr 



(1) 



where rjj is the black hole event horizon, and (3 is a con- 
stant for a given value of the black hole spin parameter a. 
It is defined by 



= 



nsco 



1 . 



(2) 



Here the radius of the innermost stable circular orbit 
(ISCO), nsco; is the same as r; n in Artemova et al. l|1996fl . 
Both rn and nsco depend on a as shown in Fig. ^ 

This potential reduces to the usual Paczyhski-Wiita 
potential (Paczynski & Wiita ITI)80|I when the black hole 
spin parameter a is zero. It allows us to mimic important 
general relativistic effects, in particular to reproduce the 
existence of a last stable circular orbit. 

The local Keplerian velocity for the Artemova- 
Bjornsson-Novikov potential is 



VKepler(r) 



lGM(r)rP-i 
(r - m)P 



(3) 



with M(r) being the mass enclosed by the sphere of ra- 
dius r and (3 being given by Eq. This becomes the 
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Table 2. Some results for the torus evolution models. The dynamical simulations were performed over a time interval 
of Ai ca i. All quantities refer to the conditions found at the end of the simulations. is the mass of the gas of the 
disk, T max is the maximum gas temperature in energy units, L„ c denotes the electron neutrino luminosity near the 
end of the simulation, L Pc the corresponding electron antineutrino luminosity, and L„ x the luminosity of heavy-lepton 
neutrinos (summed for v^, D^, v T or D T , which are all treated equally). The sum of all neutrino luminosities is given 
by L u , and the mean energies of the different neutrino types by (e Ue ), (ep c ) and (e^). 
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Keplerian velocity of the Paczyhski-Wiita potential when 
the black hole spin parameter a is zero and thus /3 = 2. 

The radius of the vacuum sphere mimicking the black 
hole is chosen to be the arithmetic mean of event horizon 
and innermost stable circular orbit. The outflow bound- 
ary condition allows mass and angular momentum (and 
energy) to leave the grid. The mass and the spin of the 
black hole are accordingly updated during the simulations. 
The corresponding changes of the black hole spin are, how- 
ever, minor (cf. Setiawan et al. 120040 because the accreted 
mass is small compared to the initial mass of the black hole 
in all cases. 

So the only GR effect being modeled is the reduction in 
radius of the ISCO. Explicit frame dragging is omitted, as 
well as the gravitational lensing effect which might capture 
some of the neutrinos emitted in this region. A detailed 
discussion of the lensing effect can be found in Sect. 7.2 of 
Ruffert & Janka lTMSjl , where we find a moderate decrease 
of the annihilation energy by 10%-30%. 



2.2. Viscosity 



In our simulations we compare the evolution of the disk 
with and without shear viscosity. The viscosity source 
terms for energy, momentum components, and entropy 
equations can be expressed, respectively, as follows 
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(5) 



(6) 



where r/ p is the dynamical shear viscosity coefficient pa- 
rameterized by the simple expression: 



r] p = ap(? s /Vt K 



(7) 



Here a is the standard dimensionless disk-viscosity pa- 
rameter introduced by Shakura and Sunyaev (19731, Q,k 
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is the Keplerian angular velocity, and c s = (TP/p) 1 ^ 2 is 
the adiabatic sound speed with T being the adiabatic in- 
dex. The expressions of Eqs.0] and were added as source 
terms to the energy and momentum evolution equations 
of our PPM scheme. 

We use the subscript p for the shear viscosity that 
is assumed to result from physical processes in the disk 
(e.g. turbulence or magnetohydrodynamic processes, see 
Balbus & Hawley [TUM Duez et al. I2UU1 Thompson et 
al. E005) to distinguish it from the numerical viscos- 
ity which was discussed in Ruffert & Janka ( 2001) and 
Ruffert & Janka I1999|) . The approach used in this pa- 
per is similar to the disk formalism developed by Chen 
et al. lfH?95). Abramowicz et al. (1995), Lee & Ramirez- 
Ruiz l|2002[) . and McKinney & Gammie (|2"UU2"|) : see also 
Hawley & Krolik lj2"UUU|) . Narayan et al. (HEB, and Kohri 
k Mineshige (|2UU12j) . 

The resolution dependent viscous dissipation due to 
numerical viscosity can be estimated to correspond to an 
a value of approximately 0.01 for our code and the cho- 
sen grid zoning (see Ruffert & Janka |2"T)0 11 and Janka et 
al. lTMH) . 

2.3. Thermodynamics 

The thermodynamics of the gas is described by using the 
equation of state (EoS) of Lattimer & Swesty (|l991). 
This allows us to follow the neutrino production and es- 
cape by a leakage and trapping scheme as detailed in 
previous papers (Ruffert et al. IT$9l)l ITM71 Ruffert & 
Janka 119991 12001fl in order to calculate the energy and 
lepton number changes by neutrino losses. This approach 
is used in the present paper, because our main aim is the 
investigation of the relevance for gamma-ray burst sce- 
narios, in particular for those where the neutrino emission 
has been suggested to provide the energy for the relativis- 
tic gamma-ray burst fireball via neutrino-antineutrino an- 
nihilation. Also the amount of mass ejection during the 
dynamical interaction and the properties of the ejected 
matter depend on the EoS. 

The density is bounded from below in our simulations 
by the assumption of an environmental medium of about 
10 s g cm~ 3 (determined by the minimum density of the 
employed equation of state table) and is more than two 
orders of magnitude below the average densities inside our 
simulated tori. 

In addition to the standard equations of continuity 
and conservation of momentum and energy, we also evolve 
an entropy equation (more information can be found in 
Ruffert & Janka H2001J1 . The terms due to viscosity are 
given in Sect. 12.21 This allows us to obtain the temperature 
from both the energy and the entropy equation, allowing 
a cross-check. For the simulations of the torus evolution 
shown in this paper, the two temperatures are in fairly 
good overall agreement, so this double-check of the tem- 
perature is redundant in the context for the simulations 
presented in this paper. Local differences between the two 



temperatures remain local and are therefore of no major 
importance for the overall hydrodynamical results. 

2.4. Models and initial conditions 

Tabled contains a list of computed models with their spe- 
cific parameters, including the values of a. We follow the 
non-stationary evolution of the toroidal disk around the 
central black hole on much longer time scale than the dy- 
namical time scale of the merger of black hole/neutron 
star systems, considering the viscous transport of angular 
momentum which depends on the assumed shear viscosity 
of the disk which we parametrise in the models. Model 
"rOO" defines a reference case without black hole rotation 
and without disk viscosity. The "al" models include disk 
viscosity but no black hole spin, models "ro" and "ri" 
consider rotating black holes but no disk viscosity, and 
models "ar" take into account non-zero values for the disk 
viscosity and BH spin. The names of the models also carry 
information about the number of zones of each grid level 
(32 or 64). 

The data for our initial setups, which consist of a 
black hole and an axisymmetric, toroidal accretion disk 
with chosen masses, were taken from one of the final 
configurations obtained in the simulations described by 
Janka et al. l(T§551 see also Ruffert et al. iflMfl IH)9T|l and 
Ruffert & Janka : 2(l(TTt for more information). Density, 
temperature, and neutron-to-proton ratio of the matter, 
which is assumed to be in nuclear statistical equilibrium, 
were azimuthally averaged. Because we varied the disk 
mass by simply scaling the density, and moreover switched 
from the Newtonian potential of the merger runs to the 
pseudo-Newtonian gravity in this work, the tori arc ini- 
tially slightly out of rotational equilibrium. Despite the 
artificial adjustment which happens during the first mil- 
liseconds of our simulations, we do not consider this as 
a major problem because the compact object merger as 
well as a subsequent collapse of the merger remnant to a 
black hole are violent processes, followed naturally by a 
relaxation phase of the mass left in a torus. 

The simulations of the present work started with a 
black hole mass of 4.017 M©. This is the mass produced 
by the final configurations of the models described above. 
We will use the model with a torus mass of 0.0478 M Q , 
as "reference model" . The mass was then also changed to 
values between 0.0120 M© and 0.1912 M© (see Table HJ 
by simply scaling the density distribution of the reference 
configuration with an appropriate factor. The initial con- 
ditions of temperature, density, and electron fraction (i.e., 
proton-to-baryon ratio) for the reference case are shown 
in Fig. □ 

The initial black hole spin parameter was varied be- 
tween values of 0, 0.3, 0.6 or 0.8 for 32-resolution models 
and between or 0.6 for 64-resolution models, respectively. 
We did not choose a closer to unity because our code is 
not fully relativistic. When a is near unity, the disk moves 
in so close to the event horizon that general relativistic 
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Fig. 2. Initial configuration at the start of the simulations for Model r00-64. The density is displayed in the orbital 
plane (left) and in two orthogonal planes perpendicular to the equator (right). It is given in gcm~ 3 with contours 
spaced logarithmically in steps of 0.5 dex. The arrows in the density plots indicate the velocity field. The temperature 
and electron fraction (Y e ) are displayed in the orbital plane. The temperature is measured in MeV, its contours are 
labelled with the corresponding values. The contours of electron fraction are spaced linearly in steps of 0.02. 



effects are extremely strong and the fluid velocities ap- 
proach the speed of light, thus making the approximative 
description of relativistic gravity and the nonrelativistic 
fluid dynamics of our code questionable. 

3. Results 

In this section we will discuss only the models calculated 
with a resolution of 64 zones on each level of the nested 
grid, so the specification "-64" in the model name, e.g. rOO- 
64, will be mostly omitted. Comparing runs with 32 zones 
and 64 zones at the same time, one finds only relatively 
small differences of the model evolution and global disk 
properties, of which some are given in Table The 3D 
simulations with 32 zones, however, could be performed 
for a time interval of 70 ms, whereas computer resources 



did not allow us to follow the evolution of the better re- 
solved models for more than 40 ms. In Table [21 we also 
list the corresponding results of the 32-zone runs, which 
include cases that we studied with high resolution, too. 

3.1. Dynamical evolution of reference Model rOO 

We begin with discussing the evolution of the disk in 
Model rOO, which will serve as reference model in this pa- 
per. Figure |2 shows the initial distributions of density, ve- 
locity, temperature, and electron fraction for this model. 
The obvious axisymmetric construction of these distribu- 
tions was described above in the section on initial condi- 
tions. Figure |3 shows these distributions after 11.2 ms. At 
this time the disk has reached a quasi-equilibrium state, 
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Fig. 3. Density distribution for Models r00-64 in the orbital plane (left) and perpendicular to it (right) at about 
11 ms after the start of the simulations. It is given in gem -3 with contours spaced logarithmically in steps of 0.5 dex. 
The arrows in the density plots indicate the velocity field. The temperature and electron fraction (Y e ) are displayed in 
the orbital plane. The temperature is measured in MeV, its contours are labelled with the corresponding values. The 
contours of electron fraction are spaced linearly in steps of 0.02. 



and its properties in Model rOO are the following: Within 
100 km of the black hole the densities are a few times 
10 10 gcm -3 , the temperatures are around 2 MeV, and the 
electron fraction is 0.04-0.06. 

Given our nested Cartesian grid, one cannot expect a 
precisely axisymmetric distribution to retain its symmetry 
when evolved with an explicit code. Some degree of numer- 
ical noise will always creep into the simulation and break 
the symmetry. This is expected to be on the scale of a cou- 
ple of zones, so on a small scale as compared to the radius 
of the disk. On the other hand, the nested grid structure, 
specifically the boundaries of the grids and the transition 
to a lower resolution outside of a central cubic volume, 
will imprint a distinct 4-fold symmetry, which mostly re- 
mains small in amplitude. So, although it is known that 



disks also suffer physical instabilities, e.g. Papaloizou & 
Pringle (|1984[) . Kleiber & Glatzel (1999), our simulations 
do not permit to connect the observed destruction of the 
initial axial symmetry with such physical instabilities. The 
numerical noise must be expected to drown all but the 
fastest growing high-order and low-order modes. 

FigureQJshows the temporal evolution of the black hole 
accretion rate and the integrated mass of the accretion 
torus. A phase of high mass accretion rate is visible from 
the outset until about t i=s 10 ms. This is a transient effect 
linked to the initial construction of the disk. This phase 
can also easily be spotted in the plot for the evolution 
of the maximum temperature in Fig El The temperature 
increases transiently. The disk settles into a quasi-steady 
state characterised by a nearly constant black hole accre- 
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Fig. 4. Mass accretion rate of the black hole and and gas mass on the grid as functions of time for Models r00-64, al3-64 
and al4-64 with increasing disk viscosity. Gas mass on the grid as functions of time for the reference Model r00-64, 
the low-mass torus Model irl-64, and the high-mass torus Models ir4-64, ir5-64, ri4-64, and ar2-64. 



tion rate at a lower value and by only slowly changing val- 
ues of the average density and temperature in the torus. 
This relaxation is more pronounced in the mass accretion 
rate of the black hole (Fig. , a fact which might indicate 
that it is the aftermath of the settling in the outer parts 
of the torus, which has now made its way to the centre. 

3.2. Influence of viscosity and varied disk mass 

The distributions of density, velocity, and temperature for 
Models al3 and al4 are shown in Fig. These have the 
same parameters as Model rOO except for an increased vis- 
cosity (see Table QJ. While both models show, in Fig. El 
a clear increase in disk temperature with increasing vis- 
cosity, the spatial distribution is different: Model al3 with 
the small value of the viscosity shows high temperatures 
between 4 MeV and 8 MeV only within a central region 
of radius 100 km (Fig. [5}, whereas Model al4 has 4 MeV 
and higher in the whole extended disk. Since the pressure 
increases with temperature, this heating by the increased 



viscosity expands the disk and reduces its density (Fig. [5J . 
Assuming otherwise equal velocities, a reduced density im- 
plies a reduced mass accretion rate. However, the large- 
scale effect of viscosity is to redistribute angular momen- 
tum from the centre of the disk out to the periphery, thus 
allowing mass to drift inward and be accreted by the black 
hole. The interplay of these two effects, the increase of the 
accretion due to viscosity on large scales and the decrease 
of the accretion rate due to a local decrease of the den- 
sity in the innermost parts of the disk, can be observed in 
Fig.0]for Models al3 and al4. In Model al3 the modest vis- 
cosity affects only the inner, most rapidly orbiting regions 
of the disk where shear effects are largest. The associated 
local central heating raises the thermal pressure and thus 
conspires to reduce the mass accretion rate. When the vis- 
cosity is further increased as done in Model al4, the global 
viscosity effect dominates, and the mass accretion rate in- 
creases beyond the values of Model rOO where no viscosity 
was added. As a consequence, the disk mass drops more 
rapidly than in the other models (Fig. |3J, leading to a 
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Fig. 5. Density and temperature distribution for Models al3-64 (left) and al4-64 (right) with disk viscosities of a = 0.01 
and a = 0.1, respectively. The density is given in gem -3 with contours spaced logarithmically in steps of 0.5 dex. The 
arrows in the density plots indicate the velocity field. The temperature is measured in McV, its contours arc labelled 
with the corresponding values. 
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Fig. 6. Maximum values of gas density (left) and temperature (right) on the grid as functions of time for Models r00-64, 
al3-64 and al4-64 with increasing disk viscosity. 
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Fig. 7. Density and temperature distribution for Models ro2-64 (left) and arl-64 (right) for a corotating black hole, 
with disk viscosities of a = 0.0, and a = 0.1, respectively. The density is given in gem -3 with contours spaced 
logarithmically in steps of 0.5 dex. The arrows in the density plots indicate the velocity field. The temperature is 
measured in MeV, its contours are labelled with the corresponding values. 



clearly lower maximum density in the disk (Fig. . The 
maximum disk temperature, on the other hand, is clearly 
largest in Model al4 (Fig.[BJ. 

The model pair ir4 and ir5, which also only differs in 
the amount of added viscosity, shows the same relationship 
with respect to the mass accretion rate as Models rOO and 
al3: In Fig. 0] one can see a smaller mass accretion for the 
model with higher viscosity (ir5) . The reasoning is similar 
to what has been described above: The temperatures are 
higher for Model ir5, but the densities distinctly lower. It 
seems that for tori with higher mass, which is the case for 
Models ir4 and ir5 relative to Models rOO, al3, and al4, the 
onset of the full effect of the increase of the mass accretion 
rate due to global angular momentum redistribution is 
offset to even higher values of the viscosity parameter. 
Increasing the disk mass by a factor of 4 does not produce 



any other major changes in flow dynamics. However the 
larger densities do have a large effect on the emission of 
neutrinos as will be described further below. 

3.3. Influence of black hole rotation and of a more 
compact gravitational potential 

We mimic some effects of relativistic gravity by using a 
Paczyhski-Wiita (1980) potential and its extension to ro- 
tating black holes by Artemova et al. (|I996). Varying the 
spin parameter between negative (counter-rotating black 
hole) and positive (black hole corotating with disk) values, 
the radius of the horizon changes and the ISCO is shifted 
to smaller values compared to the case of a black hole 
without rotation (Fig. [TJ. Looking at the central circle in 
Fig.|7|as compared to, e.g., Model rOO, one sees this slight 
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Fig. 8. Left panels: The dots represent the azimuthal velocities v ip (d) (normalised to the speed of light) of all zones 
in the equatorial plane of the reference Model r00-64, and the high disk-viscosity Model arl-64 with corotating black 
hole, at the times after the start of the simulation indicated in the plots. The thin solid line gives the average value 
of the azimuthal velocities and the bold solid line the local Keplerian velocity WKepier(d) of the Artemova et al. (1996) 
potential as a function of the equatorial radius d (Eq. |2J. The inner vacuum boundary of the computational grid, is 
located at the arithmetic mean of the event horizon and innermost stable circular orbit (ISCO) of the central black 
hole. Right panels: The dots give the radial velocities v r (d) (normalised to the speed of light) as a function of the 
equatorial radius d for all grid zones in the equatorial plane of the reference Model r00-64 and the high disk-viscosity 
Model al4-64, at the same times. The solid lines represent the mean values of all zones within binning intervals of 
3 km. 



shrinking of the inner grid boundary (at the arithmetic 
mean of rjj and nsco) when the black hole rotates in the 
prograde direction. Models ro2 and ro5 can be directly 
compared to rOO, the only difference being the presence of 
rotation. The same applies to the pair al4 and arl. 

No major difference can be seen in the density or tem- 
perature distributions, except at the very centre close to 
the black hole. The temperatures for the rotating black 
holes are about 2 MeV higher than in the non-rotating 
cases The counter-rotating Model ro5 retains the low tem- 
peratures of the reference Model rOO, i.e. around 2 MeV. 
In case of the viscous tori with corotating black hole the 



higher temperatures in the close vicinity of the black hole 
correlate with a distinct narrowing of the accretion torus 
in the vertical (z) direction. 

The combined increase of central temperatures 
brought about by viscosity and rotation has a clear im- 
pact on the velocities in the direct vicinity of the black 
hole. Fig. [S] shows separately the azimuthal and radial ve- 
locities in the equatorial plane versus the distance from 
the centre of the black hole. Each zone is represented by a 
dot, and the spread of the points at a given radius reflects 
the deviations from rotational symmetry. For an axially 
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Fig. 9. Density and temperature distribution for high-mass torus models with nonrotating black hole and without 
disk viscosity (ir4-64; left) and with corotating black hole and disk viscosity a — 0.1 (ar2-64; right) at about 20 ms 
after the start of the simulations. The arrows indicate the velocity field. The density is given in gem -3 with contours 
spaced logarithmically in steps of 0.5 dex. The temperature is measured in MeV, its contours are labelled with the 
corresponding values. 



symmetric configuration all dots at a specific radius would 
cluster on top of each other. 

Comparing the radial velocities of Models al4 and rOO 
(top and bottom right panels in Fig. |SJ , one sees that the 
large viscosity produces the inward facing accretion flow, 
i.e. negative velocities, at larger radii, inside of 60 km for 
Model al4 as opposed to 40 km for Model rOO. This re- 
flects the canonical viscosity effect of globally shifting an- 
gular momentum outward allowing mass to fall inward. 
The added effect in Model arl of narrowing the potential 
well with a rotating black hole does not change this sit- 
uation significantly. Note that the outer parts of the tori 
show slow expansion at the displayed times, indicating 
ongoing relaxation as well as the consequence of outward 
transport of angular momentum. 



The azimuthal velocities in Model rOO closely follow 
the Keplerian values for the Paczyhski-Wiita potential, see 
the top left panel in Fig. [SJ Although the added viscosity 
in Model al4 (bottom panel) increases the temperature in 
the disk, the pressure change due to this increase is not 
sufficient to appreciably change the velocity distribution. 
Only when black hole rotation is added, with its additional 
temperature increase and corresponding pressure increase 
in the vicinity of the black hole, does the matter in the disk 
become more strongly pressure supported within a radius 
of 50 km. The azimuthal velocities are significantly below 
the Keplerian values in this region, as given by Eq. 
(with the radius r replaced by the equatorial distance d) 
for the Artemova-Bjornsson-Novikov potential. 
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In all four pairs of models, r00/ro2, ri4/ir4, al4/arl, 
and ir5/ar2 (see Figs. ED for snapshots of the most massive 
torus Models ir4 and ar2 at t k 20 ms) the first mem- 
ber of the pair is the non-rotating case, while the second 
contains a corotating black hole. All other parameters are 
kept equal. Examining Tabled one notices that for each 
case the rotating model has a higher disk mass than the 
corresponding non-rotating model at the end of the sim- 
ulated evolution (cf. also Fig. ^J. This is a strong indi- 
cation that a corotating black hole accretes matter from 
the initial tori less quickly, because the angular momen- 
tum at the ISCO is lower and a rotating black hole thus 
allows matter with lower angular momentum to remain 
on orbits. Moreover, the higher temperatures close to the 
black hole provide pressure support and additionally sta- 
bilise the torus. A confirmation is found by looking at 
the counter-rotating case ro5, which can be compared to 
Models rOO and ro2. Indeed, the disk of Model ro5 has an 
even smaller mass than that of Model rOO. 

3.4. Neutrino emission 

The neutrino luminosities and cumulative energy emit- 
ted in neutrinos as functions of time for a sample of 64- 
resolution models are shown in Fig. Note that electron 
neutrinos, v e , and electron antincutrinos, v ei dominate the 
emission by a large factor, because at the density and 
temperature conditions in the accretion tori their produc- 
tion rate by charged-current absorption of electrons and 
positrons on free protons and neutrons, respectively, 

e~ + p — > n + v e , (8) 
e + + n — ► p + D e , (9) 

is much higher than the production rate of neutrino- 
antineutrino pairs of all flavors by thermal neutrino 
processes and nucleon-nucleon bremsstrahlung (see also 
Ruffert & JankaHnHHI Rosswog & Liebendorfer EMU ■ 

The v e luminosities are initially much larger than the 
v e luminosities and at the end of the calculated evolu- 
tion still typically 2-3 times higher (see Table |2J). The 
tori consist of decompressed neutron star matter, which 
inside of the neutron star was initially in a very neutron- 
rich (Y e ^0.1) state of neutrinolcss beta equilibrium (the 
initial setup chosen for our simulations was adapted to 
neutron star merger results and therefore still reflects this 
origin of the torus matter, cf. Fig. Chemical or ki- 
netic equilibrium at densities much below nuclear matter 
corresponds to larger values of Y e and a higher proton 
abundance. At decompression, matter tries to approach 
this new equilibrium state, in which the production of v e 
and v e are in balance. During this "protonization" electron 
antincutrinos are emitted in larger numbers than v e . The 
consequences of this process are visible in two ways. First, 
comparing the Y e distribution in Fig. ^2 after 11ms and 
20 ms of evolution with the initial state shown in Fig. 
one sees significantly higher values of Y ei i.e. more pro- 
tons. Secondly, in Fig. ^1 the comparison of the local 



energy emission rates per unit area or volume, of u e , D e , 
and heavy- lepton neutrinos plus antincutrinos ( u v x " f° r 
the sum of muon and tau neutrinos, Vp, v^, v T , and v T ), 
respectively, shows that the v e emission is significantly 
higher than that of v e . The protonization is faster and 
leads to higher values of Y e at the end of our simulations 
in case of hotter and less dense tori; this trend is very simi- 
lar for larger viscosity and black hole rotation (Models ro2, 
arl, ir4 and ar2 in Fig. 1111) . All tori are at most marginally 
opaque to neutrinos in their denser parts and the diffusion 
times short in most regions. The Y e evolution is therefore 
governed by the capture rates directly. Since positron cap- 
tures in reaction increase extremely steeply with rising 
temperature, the highest values of Y e develop in regions 
of highest temperatures. 

The mean energies of the radiated neutrinos ((ej/ e ), 
(ej/ e ) , and (e„ x ) ) are defined as the ratio of integral energy 
loss (rate) to total number loss (rate) in neutrinos of a 
given kind (see Ruffert et al. [1996 1997}. They correlate 
with the average and maximum torus temperatures as visi- 
ble from the parallel trends of T max and the mean neutrino 
energies in Table |21 Electron neutrinos and antineutrinos 
are radiated away with mean energies of 10-20 MeV and 
are similarly energetic, because electron degeneracy plays 
a minor role and the cross sections of both production 
reactions scale equally with particle energies; moreover, 
reabsorption of neutrinos is of minor importance because 
the tori are essentially transparent to neutrinos. Heavy- 
lepton neutrinos, on the other hand, tend to have a slightly 
lower mean energy except in cases with the largest torus 
masses (Models ri4, ir4, ir5, and ar2), which also tend to be 
among the hottest models. Only in the two cases with the 
highest torus temperatures and shear viscosity parameter 
a = 0.1 (Models ir5 and ar2) do muon and tau neutrinos 
yield a non-negligible contribution to the total energy loss 
in neutrinos. 

The results summarised in Tables [5] and [3] and Fig. 1101 
show that the neutrino luminosities rise steeply with torus 
mass. In simulations without shear viscosity (a — 0) we 
find a roughly quadratic increase, L v cx M 2 (compare 
Model r00-64 with ir4-64 and ro2-64 with ri4-64), whereas 
with shear viscosity a = 0.1 a certain saturation can 
be observed and the luminosities increase only slightly 
steeper than linearly, L v cx M 1+ ^ for £ > (Models al4- 
64 vs. ir5-64 and arl-64 vs. ar2-64). Shear viscosity has 
the most dramatic influence on the neutrino emission: the 
luminosities grow by a factor of 30-40 when a — 0.1 in- 
stead of a = (Models r00-64 vs. al4-64; ir4-64 vs. ir5- 
64). In contrast, a black hole with a spin parameter up to 
a ~ 0.6 in corotation with the torus has a much weaker 
influence on the neutrino emission. While for calculations 
without shear viscosity (a — 0) the neutrino luminosities 
increase at most by a factor of 2-3 (compare Model r00-64 
with ro2-64 and Model ri4-64 with ir4-64 in Table □ and 
Fig. 110(1 , an additional increase of the neutrino emission 
is hardly visible or absent for cases with a = 0.1 (com- 
pare Models al4-64 with arl-64 and Models ir5-64 with 
ar2-64). Black hole rotation with more extreme spin pa- 
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Fig. 10. Total neutrino luminosities (left) and cumulative energy radiated in neutrinos (right) as functions of time for 
the reference Model r00-64, the low-mass torus Model irl-64, and the high-mass torus Models ir4-64, ir5-64, ri4-64, 
and ar2-64. 



rameter, however, leads to a large boost of the neutrino 
emission as can be seen from Models ro2-32 with a = 0.6 
and ro3-32 with a = 0.8 in Table [5] Despite of nearly 
equal torus masses, the latter model has a 9 times higher 
neutrino luminosity. 

The series of models with reduced resolution exhibits 
the tendency of slightly lower torus temperatures com- 
pared to the corresponding simulations with finer grid 
zoning at the same time of the evolution. The neutrino 
luminosities and mean energies, however, show no clear 
trend and can be higher or lower, depending mostly on 
whether the torus has a larger or smaller mass than in 
the better resolved calculations. This sensitivity together 
with the fact that the less massive tori have lost most of 
their mass to the black hole after an evolution of 70 ms 
(Table disfavours the late-time models with lower res- 
olution for a quantitative comparison of the influence of 
viscosity and black hole rotation. The neutrino emission 
during these late stages has dropped to a low level in most 
cases and is therefore mostly dependent on the torus mass 
which is left around the black hole (see Table [2J . The 
general trends associated with viscosity, torus mass, and 
black hole rotation, however, can also be seen in the 32- 
resolution models. 

The upper left panels of Fig. ^5] and right panel 
of Fig. 1141 provide information about the instantaneous 



conversion efficiency of rest-mass energy to neutrino- 
antineutrino energy, 



for a selection of 64-resolution models as function of time. 
It is obvious that the values are extremely low (q v <C 
1%) in case of zero physical shear viscosity (a — 0) and 
become of order 1% only for large torus masses (Ma > 
0.1 M Q as in case of Model ri4-64, cf. TableEJ). For a = 0.1 
the efficiencies reach several percent (e.g. Model arl) and 
Models ar2 and ir5 in Fig. I13|) . It is interesting that q v 
for the low-mass torus Model arl reaches its maximum 
at early times and then decreases towards the end of our 
simulations, whereas the values are stable or even increase 
over the simulated period of evolution in Models ir5 and 
ar2. 

3.5. Neutrino-antineutrino annihilation 

Energy deposition in the vicinity of the black-hole torus 
system by the annihilation of neutrinos (v) and antineu- 
trinos (u) of all flavors, which are radiated away from the 
accretion torus, to electron-positron pairs in the process 

v + v — ► e + e~ (11) 

is considered as a way to produce a highly relativistic pair- 
plasma, if the density of baryonic matter in the deposition 
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Fig. 11. Electron fraction in the orbital plane for Models ro2-64, arl-64, ir4-64 and ar2-64, in the orbital plane at 
about 11ms after the start of the simulations. The contours are spaced linearly in steps of 0.02. 



region is sufficiently low. Provided the entrainment of ad- 
ditional baryons during the subsequent expansion can be 
prevented, the pressure by radiation and electron-positron 
pairs can then accelerate the plasma to ultrarelativistic 
velocities. Recent general relativistic hydrodynamic simu- 
lations have demonstrated that this indeed happens and 
that the ultrarelativistic outflow becomes collimated to 
semi-opening angles between ~ 5° and ^12° (Aloy et 

ai.imi. 

Because we calculate energy and lepton number loss 
in our 3D hydrodynamic simulations of the torus evolu- 
tion by a neutrino trapping scheme, but do not solve the 
neutrino-transport problem, the corresponding energy de- 
position is evaluated for our hydrodynamic models in a 
post-processing step. The numerical procedure for doing 
this is rather CPU-time consuming and explained in detail 
in Ruffert & Janka 1(13951 ISUUTj) and Ruffert et al. (fTW7|) . 
The energy deposition rate per unit volume, Q^p, at a 



point r is calculated as the sum of the contributions from 
neutrinos and antineutrinos emitted from all regions of 
the accretion torus, thus taking into account the specific 
emission geometry and local variations of the neutrino loss 
rates as visible in the emission maps of Figs. ll2| Neutrinos 
and antineutrinos radiated from the hot accretion torus 
interact with each other in the surroundings with a finite 
probability, which depends on the number densities and 
the energies of these neutrinos and on the angle between 
the directions of neutrino and antineutrino propagation 
(see, e.g., Goodman et al. HMTl Ruffert et al. HMTjl . The 
total energy deposition rate, 



E, 



C L v Li, 



(12) 



therefore increases with the product of neutrino and an- 
tineutrino luminosities and with the spectrally averaged 
neutrino energies, times a factor C that accounts for the 
dependence on the angular distribution of the neutrinos. 
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Fig. 12. Neutrino energy loss rates per unit area in the orbital plane for Model r00-64 at 20 ms after the start of 
the simulations. The plotted values show the logarithm of the rates in ergcm _2 s _1 , obtained by integration of the 
local energy loss rates per unit of volume from z = to infinity. The top left panel gives the results for electron 
neutrinos, the top right panel for electron antineutrinos, the lower left panel for the sum of muon and tau neutrinos 
and antineutrinos, and the lower right panel the total values for neutrinos and antineutrinos of all flavors. The contours 
are spaced in steps of 0.5 dex, bold lines are labelled with their corresponding values. The grey shading emphasises 
the emission levels, dark grey representing the largest energy loss by neutrino emission. 



The quantities (e„) and (e 2 ) in Eq. 1112[| denote the mean 
energies and mean squared energies of v and v. The factor 
C contains the weak interaction coefficients and terms that 
depend on the geometry of the neutrino-emitting torus 
region. This factor is a function of the emission geome- 
try and drops steeply with the distance from the neutrino 
source. 

Large collision angles between neutrinos and antineu- 
trinos are found only close to the source and therefore an- 
nihilation is favoured near the hottest parts of the torus, 
close to the equatorial plane. There, however, the gas den- 
sities are still high and the energy deposition by vv anni- 
hilation is compensated by very rapid cooling due to the 



charged-current reactions of Eqs. JHJ) and Q 1 . Neutrino- 
antineutrino annihilation is therefore not efficient in es- 
tablishing net energy transfer to regions with large baryon 
densities. On the other hand, the annihilation rate is also 
higher above the poles of the black hole, where neutrinos 
emitted from the torus have a high probability to collide 
with large angles in low-density environment. This can be 
seen in the four panels of Fig. which show the an- 
nihilation rate in the d-z-planes perpendicular to the or- 



1 The energy loss associated with the reemission of neutri- 
nos by charged-current reactions is not taken into account in 



Fig. Ha 
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Fig. 13. Conversion efficiency of rest-mass energy to neutrinos, q v = L„/(Mc 2 ) (top left), conversion efficiency of 
neutrino energy to electron-positron pairs by neutrino-antineutrino annihilation, q vi/ = E vi >/L u (top right), integral 
rate of energy deposition by neutrino-antineutrino annihilation around the accretion torus, E^c, (bottom left), and 
cumulative energy deposition by vv annihilation, E V y (bottom right), as functions of time for the same set of models 
shown in Figs. 01 and 1101 
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Fig. 14. Integral rate of energy deposition by neutrino-antineutrino annihilation around the accretion torus, E v „ 
(left) and conversion efficiency of neutrino energy to electron-positron pairs by neutrino-antineutrino annihilation 
g„ p = E V vjL v (right), as functions of time for selected 64-resolution models with the same torus mass but different 
torus viscosity. The models are the same as in Figs. 0] and EJ 



Table 3. Torus mass Md, black hole mass accretion 
rate Md, estimated accretion time scale of the torus, 
i acc = Md/Md, total neutrino luminosity L v , integral rate 
of energy deposition by neutrino-antineutrino annihilation 
around the accretion torus, E u p, and total energy deposi- 
tion by vv annihilation in the time interval tfi n , E v „, for 
all 64-resolution models. All quantities are given at time 
£fi n after the start of the simulation. 



Model 



£fin 

ms 
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10~ 2 M o 
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bital plane for Models rOO-64, ro2-64, ro5-64, and ar2-64 
at 20 ms after the start of the simulations. 



The corresponding integral rate of energy deposition 
by neutrino-antineutrino annihilation around the accre- 
tion torus, E v „, the cumulative energy deposition by vv 
annihilation, E V y, and the neutrino conversion efficiency, 

q u y EE Evv/Ly , (13) 

as functions of time are plotted for some of our computed 
models in Figs. El an d El Values of E v p and E u p for 
all 64-resolution models at the end of the simulations are 
listed in Table 

The contributions from muon and tau neutrinos and 
antineutrinos can safely be ignored in this evaluation, be- 
cause their luminosities are much lower (see Table |2J and 
their annihilation cross sections are also smaller by a fac- 
tor of about 5 compared to those of electron neutrinos 
and antineutrinos. The geometry of the neutrino-emitting 
torus region and therefore the factor C does not change 
strongly during the quasi-stationary phase of the torus 
evolution and therefore C can be calibrated by evaluating 
the vv annihilation in detail at one or more representative 
times. Thus Eq. p2|) with constant C allows for a good 
approximation of the time evolution of E u p (as shown in 
Figs. ITS! and ITU) by using the neutrino luminosities and 
mean neutrino energies as determined from the torus sim- 
ulations, without the need to perform the CPU-time de- 
manding post-processing for the ^-annihilation maps of 
Fig. 1 151 many times for our sample of runs. 
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Fig. 15. Maps of the local energy deposition rates (in ergcm _3 s _1 ) by vv annihilation to e + e" pairs in the sur- 
roundings of the accretion torus for Models r00-64, ro2-64, ro5-64, and ar2-64 at a time 20 ms after the start of the 
simulation. The solid contour lines in a plane perpendicular to the equatorial plane represent values averaged over 
azimuthal angles around the z-axis. The contours are logarithmically spaced in steps of 0.5 dex and the grey shading 
emphasises the levels with dark grey meaning high energy deposition rate. The energy deposition rate was evaluated 
only in that region around the torus where the baryon mass density is below 10 11 gcm~ 3 . The white octagonal area 
around the centre with a semidiameter of one Schwarzschild radius indicates the presence of the central black hole. The 
dotted contour lines represent levels of constant values of the azimuthally averaged mass density, also logarithmically 
spaced with intervals of 0.5 dex; the bold dotted line corresponds to p = 10 10 gcm~ 3 . 



Only our most massive tori possess initially neutrino- 
opaque "cores" with densities above ~ 10 11 gcm _1 , 
whereas in the other cases such regions are absent. When 
the optical depths for neutrinos are larger than unity and 
ncutrinospheres exist, these typically have ellipsoidal or 
egg-shaped cross sections in the d-z plane. Also in our 
massive tori the opaque cores disappear and the neutrino 
optical depths fall below unity everywhere after the tori 
have lost some mass to the black hole and/or thermal in- 
flation by viscous heating has decreased the density. All 



of the tori shown in Fig. 1151 are therefore transparent 
to neutrinos. A comparison of the annihilation maps of 
Models rOO, ro2, and ro5 in Fig. 1151 shows the influence of 
the black hole rotation. The values of the energy deposi- 
tion rates differ and scale roughly with the neutrino lumi- 
nosity as expected from Eq. i|12|) (cf. also Table|3J), but the 
morphology of the contour lines which envelope the black- 
hole torus system is very similar. In the counter-rotating 
case (Model ro5) the two central maxima are more pro- 
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nounced and represent a toroidal volume which is more 
compact than in the other cases. 

Only our most massive tori (Models ri4, ir4, ir5, and 
ar2) and the viscous intermediate-mass Models al4 and arl 
reveal vv annihilation rates E vij which become larger than 
10 50 ergs _1 at least for some period of the simulated evo- 
lution (see Figs. 1131 and !14l the non- viscous massive torus 
Model ir4, however, only marginally reaches this value). 
The numbers listed in Table |3| clearly show the correla- 
tion with the neutrino luminosities. They also show that 
the black hole rotation (up to a rotation parameter of 
a = 0.6) has hardly any influence on the total energy E u „ 
deposited by vv annihilation (compare the model pairs 
rOO and ro2, al4 and arl, and ir5 and ar2). In contrast, 
the torus mass makes a big effect for viscous (model pairs 
al4 and ir5, arl and ar2) as well as nonviscous cases (model 
pairs rOO and ir4, ro2 and ri4); in case of the higher torus 
masses E vP is larger by factors of 5-40. Again shear vis- 
cosity (a = 0.1 instead of a — 0) makes the biggest differ- 
ence. For the latter, viscous models of the pairs (ri4,ar2), 
(ir4,ir5), and (r00,al4) E vV is higher by factors between 
15 and 1000. The lower right panel of Fig. 1151 shows the 
annihilation map for one of the two most extreme cases, 
Model ar2, at 20 ms after the start of the simulation when 
the integral rate of energy conversion by vv annihilation 
was E uj} — 5.8 x 10 51 ergs -1 . The energy deposition rates 
above the poles of the black hole reach 10 30 ergcm -3 s _1 at 
this time. Even at z-distances of 50 km they are still more 
than 10 29 ergcm _3 s _1 . For the most massive of our tori 
and the viscous tori with intermediate masses, the anni- 
hilation efficiencies <jw can reach several percent, peaking 
for short times at more than 5% (Tigs. H31 and IT4*|l. We 
note that the results for vv annihilation which we have 
obtained for the 32-resolution runs are in general smaller 
by about 20% compared to the better resolved models. 
We also point out that the total energy deposition by 
vv annihilation {E w ^) of Model ir5-64 at the end of our 
simulations (at 40 ms) is about 25% higher than that of 
Model ar2-64, but we estimate the remaining torus life- 
time of Model ar2-64 to be more than two times longer. 

At the end of the simulated evolution, the tori have lost 
typically about half of their mass to the black hole, and the 
mass accretion rates have decreased to 0.1-0.4 M s _1 in 
most cases. The total neutrino luminosities have declined 
from their peak values by roughly an order of magnitude 
for the tori with intermediate mass and by about a fac- 
tor of 5 for our most massive cases fFig. 1101 and Table 0). 
Correspondingly, the vv annihilation rates have dropped 
by roughly two orders of magnitude below their maxi- 
mum values in the former case (Fig. I14[) . while in the 
latter case they have decreased only by a factor of 5-10 
(Models ar2 and ir5 in Fig. ^3 which have peak values 
of E vjj ^ 10 52 ergs _1 ). This is less steep than the square 
of the total neutrino luminosities and suggests a time- 
dependent change of the emission of electron neutrinos 
relative to electron antineutrinos (note the dependence of 
E V y on the product L^Lp in Eq. I12|) . This fact is con- 
sistent with the large rise of Y e relative to the initially 



low values (e.g. visible in the comparison of Fig. 1111 with 
Fig. |5J). When Y e climbs closer to 0.5 and thus tries to 
approach a value where kinetic equilibrium of v e and v e 
production holds, the initially much stronger emission of 
v e evolves to a more balanced release of electron neutrinos 
and antineutrinos 2 . 

4. Axial low-density funnel and GRB production 

Near the poles of the black hole and along the system axis, 
the density in our simulations decreases quickly. Material 
which is swept into these polar regions during and im- 
mediately after the merging of the double neutron stars 
or neutron-star black-hole binary falls into the (newly 
formed) black hole within a free-fall time scale, because it 
is not supported by a strong negative pressure gradient or 
centrifugal forces. Thus a cylindrical funnel along the ro- 
tation axis of the black-hole torus system is rapidly "evac- 
uated" . The continuous entrainment of more gas from the 
torus into that region is hindered by the centrifugal bar- 
rier which this gas experiences due to its large angular 
momentum. 

Attempts to estimate the characteristic properties of 
ultrarelativistic outflows (e.g., its terminal Lorentz fac- 
tor or the collimation factor) on the basis of results from 
hydrodynamic merger or accretion simulations like the 
ones presented here must be taken with great caution and 
might be quantitatively meaningless. The formation of col- 
limatcd outflows in the vicinity of a torus-accretor system 
is a violent, time-dependent hydrodynamic process. The 
jet could clean its own funnel through an initially baryon- 
loaded environment as happens in the core of a collapsar, 
where large amounts of gas envelope the black-hole and 
accretion disk. On the other hand, Kelvin-Helmholtz in- 
stabilities due to the interaction of the shear flow with 
the dense torus may lead to baryon entrainment and col- 
limation. Reliable and quantitatively meaningful answers 
therefore require the modelling of the jet formation and 
evolution by hydrodynamic simulations (for more details, 
see Aloy et al. EDDS) . 

The low-density polar funnels along the rotation axis 
of the torus are unlikely to form near the poles of a hot, 
neutrino-radiating massive neutron star, if such a merger 
remnant escapes the collapse to a black hole. At the sur- 
face of such a neutrino-emitting compact object, a dense 
baryonic wind is driven by neutrino energy deposition. 
This wind produces an expanding cloud of baryonic mat- 
ter around the star, a phenomenon which is well-known 
for the hot proto-neutron stars born in supernova ex- 
plosions (Duncan et al. 113551 Woosley I1993bl Q ian & 
Woosley [EnSI Thompson et al. 127)171)1 . 

The hot merger remnant emits neutrinos and antineu- 
trinos with very high luminosities. These energetic neutri- 
nos diffuse and then stream outward to the surface where 
they must deposit energy by charged-current absorption 

2 Note that the product L v La has a maximum for L v — L& 
if L v + La is a constant. 
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on free nucleons (the inverse of the reactions Eqs. © 
and 10) in the layers of decreasing temperature outside 
of the neutrinosphere. This energy deposition drives the 
subrelativistic, baryonic wind, whose formation cannot be 
treated in existing neutron star merger simulations, be- 
cause the simulations do not include neutrino transport. 
Neutrino emission is at best described by a neutrino- 
trapping scheme, in which neutrinos arc released locally, 
but the transport of energy by the neutrinos from lay- 
ers deep inside the merger remnant to the surface is not 
followed. Because of the absence of the neutrino-driven 
wind, the densities in the polar regions of the massive 
post-merger neutron stars are largely underestimated by 
the current hydrodynamical models. 

The baryon "pollution" by the neutrino-driven wind 
makes it unlikely that such a hot neutron star merger 
remnant can be the central engine of a cosmic gamma- 
ray burst, in particular if the energy release that pow- 
ers the GRB is considered to happen in the close vicin- 
ity of the neutron star. The baryon loading in this re- 
gion will prevent the formation of a pair-plasma fireball 
or jet which can accelerate to ultrarelativistic velocities. 
Neutrino-antineutrino annihilation, for example, deposits 
energy efficiently only very close to the neutrinosphere 
and the annihilation rate drops extremely steeply with 
distance (e.g., like r~ s if the neutrino source is a spheri- 
cal or ellipsoidal body). In case of a hot neutron star the 
gas density outside of the neutrinosphere is very high — 
in contrast to the baryon-poor conditions above the poles 
of a black hole — and therefore vv annihilation competes 
with the charged-current absorption of neutrinos by nucle- 
ons (Eqs. (jHJ) and ©) and with the very efficient cooling 
due to the inverse of these reactions. In detailed wind mod- 
els it is therefore seen that vv annihilation yields only a 
minor contribution to the energy transferred to the bary- 
onic wind (see, e.g., Fig. 6 in Thompson et al. 2001), and 
this additional energy deposition triggers enhanced mass 
loss from the surface of the neutron star instead of driving 
ultrarelativistic expansion. Similar arguments also apply 
for models where the GRB outflow is considered to be 
powered by magnetohydrodynamic activity of a rapidly, 
differentially rotating hot neutron star. Also in this case 
it has to be demonstrated that the high baryon densities 
associated with the neutrino-driven wind around the neu- 
tron star do not represent a unsurmountable obstacle for 
the development of ultrarelativistic fireballs or jets. 



5. Summary 

We have presented 3D hydrodynamic simulations of the 
time-dependent evolution of black-hole accretion-torus 
systems with parameters typical of the remnants of binary 
neutron star or neutron-star black-hole mergers. The tori 
which are created from neutron stars disrupted in such 
violent events, and the conditions in the torus gas are not 
well described by steady-state assumptions and the gas 
evolution requires time-dependent modelling. 



The black hole was described by a vacuum boundary 
condition and pseudo-Newtonian gravitational potential, 
corresponding to an initial mass of about 4M . The ef- 
fects of black hole rotation were included according to 
Artemova et al. ; 1 !)!)(>). taking into account the variation 
of black hole mass and spin by the accretion of matter 
from the torus. The torus gas was treated with Newtonian 
dynamics and self-gravity, its shear viscosity was mod- 
elled by a simple a-law, and the initial conditions were as- 
sumed to be cool (temperatures of 1-2 MeV) and neutron- 
rich (proton-to-baryon ratio Y e ~ 0.02-0.1), decompressed 
neutron star matter in nuclear statistical equilibrium. The 
initial spatial distribution of the gas was modelled by us- 
ing results from (Newtonian) merger simulations. 

Assuming different initial masses of the torus (Afd.ini ~ 
0.01-0.2 Mq), black hole spins (spin parameter a = 
0, 0.3, 0.6, 0.8), and gas viscosity (a-parameter between 
and 0.1), we followed the accretion process and the as- 
sociated neutrino emission for 40 ms with high-resolution 
runs, and for 70 ms with lower resolution runs. The global 
parameters of both simulations at the same times agree 
satisfactorily well, but the better resolved simulations ex- 
hibit slightly higher gas temperatures and correspondingly 
larger neutrino emission. 

Since our initial conditions with chosen torus mass and 
black hole potential did not correspond to exact rotational 
equilibria, the first 10 ms of our simulations are charac- 
terised by a relaxation phase with high mass accretion 
rates and a dynamical reconfiguration of the torus mat- 
ter, followed by a more "quiet", quasi-stationary phase in 
which the mass accretion rate of the black hole settles to 
a lower, slowly changing value. 

The tori start to evolve quickly from their initially 
neutron-rich state to more proton-rich conditions. This 
protonization proceeds by the emission of electron an- 
tineutrinos, which dominates over the emission of elec- 
tron neutrinos and heavy-lepton neutrinos, and is faster 
for hotter tori. The gas viscosity thus has an extremely 
strong influence on the torus evolution. While simulations 
without physical viscosity (a — 0) show only compres- 
sional heating and thus stay relatively cool with Y e chang- 
ing rather slowly, shear viscosity (a = 0.01 or 0.1) leads 
to a rapid rise of the gas temperature, enhanced accretion 
of the torus regions close to the black hole, and thermally 
driven inflation of the outer parts of the torus with con- 
sequently lower densities. Therefore the mass of viscous 
tori can decrease initially more rapidly and then contin- 
ues to decline more slowly with a much reduced mass-loss 
rate to the black hole. After 40 ms of evolution the tori 
have lost about half of their initial mass to the black hole, 
which at that time accretes matter at rates between about 
0.2 Mq s" 1 and about 0.6 M & s _1 in all simulated cases ex- 
cept those with a = and the highest and lowest initial 
torus masses. 

The maximum densities in the tori have dropped to 
^ 10 11 gcm -3 after 10 ms in most cases, and stay around 
2 x 10 11 gcm~ 3 after the early relaxation phase only in 
our most massive accretion tori. The tori are therefore 
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optically thin to neutrinos or develop neutrinospheres just 
marginally. For this reason neutrinos, once created, stream 
off essentially freely without being hindered by scattering 
and reabsorption. Due to the reduced gas densities and 
relatively high temperatures in viscous tori, the protoniza- 
tion there proceeds quickly and produces Ye-values of 0.3- 
0.4 within only 20 ms, tending to establish kinetic equi- 
librium conditions for the competing v e and v e emission 
processes. As a consequence, the v e and v e luminosities 
become more similar toward the end of our simulations, 
although those of electron antineutrinos are still a factor 
of 2-3 higher. Muon and tau neutrinos and antineutrinos 
contribute to the energy loss of the accretion tori on a 
negligible level, because their creation by nucleon-nucleon 
bremsstrahlung and electron-positron pair annihilation is 
much slower than the charged-current production of v e 
and v e at the temperatures and low densities of the tori. 

The effect of the protonization on the neutrino lumi- 
nosities has interesting consequences for the rate at which 
neutrinos and antineutrinos deposit energy around the 
black-hole torus system by the annihilation to electron- 
positron pairs. While the total neutrino luminosity is 
strongly correlated with the torus mass and therefore 
tends to decay with falling gas mass, the energy deposition 
scales with the product of neutrino and antineutrino lumi- 
nosities, L v Ly (see Eq.^J electron neutrinos and antineu- 
trinos make by far the dominant contribution) . It therefore 
decreases much less steeply with time than the square of 
the total neutrino luminosities, because the electron neu- 
trino emission is initially very small but rises continuously 
relative to the one of electron antineutrinos. 

As mentioned above, shear viscosity has the largest in- 
fluence on enhancing the neutrino production. We found 
the highest luminosities and largest vv annihilation rates 
for our simulations of viscous tori (a — 0.1). Values in ex- 
cess of 10 53 ergs _1 and up to ^ 10 50 ergs _1 for L v and 
-E„p, respectively, were reached even for relatively low 
torus masses on only some 10~ 2 Mq. The maximum con- 
version efficiencies of rest-mass energy to neutrinos could 
reach 5-10%, and of neutrinos to electron-positron pairs 
3-5%. While viscosity, however, tends to decrease the 
torus lifetime (measured as the ratio of torus mass to mass 
accretion rate) , black hole rotation in the direction of the 
torus spin has the opposite effect. The higher disk tem- 
peratures close to the more compact ISCO in this case, 
however, have a noticeable effect on the neutrino emission 
only if the gas viscosity is low or if the black hole spin is 
close to maximal (a <; 0.8). 

An evaluation of our results for the growth of nonra- 
dial instabilities, the presence of which might be suggested 
by the fact that we found the axial symmetry of our ini- 
tial models destroyed in a few revolutions, is prevented 
by our use of a Cartesian grid with three nested levels of 
different resolution. This grid introduces visible numerical 
artifacts, in particular imprints a distinct 4-fold symmetry, 
and must be suspected to be also the reason for higher- 
order azimuthal modes. 



6. Conclusions and outlook 

Our models are another contribution to the ongoing efforts 
to theoretically study the viability of accretion tori around 
stellar-mass black holes as central engines of cosmic GRBs. 
A detailed, critical assessment of previous attempts to link 
the results of similar simulations to GRB properties has 
already been made in Sect. El The crucial results of our 
studies in this respect are the neutrino luminosities and 
the energy conversion by neutrino-antineutrino annihila- 
tion, the former being the direct output from the neu- 
trino trapping/leakage scheme applied in our model cal- 
culations, the latter being obtained by a post-processing 
of our results. A direct calculation of the energy release 
by ^-annihilation would require the application of a 3D 
transport description, which remains a true challenge for 
the future. 

The results of our studies show that for values of the 
torus mass, shear viscosity, and black hole rotation well 
in the range of expectations for compact binary mergers, 
the neutrino emission is sufficiently powerful and the as- 
sociated efficiencies for converting gravitational binding 
energy to neutrinos as well as neutrino-antineutrino pairs 
to a e + e _ -photon plasma, are sufficiently large to account 
for the energetics of short GRBs, in particular if colli- 
mation of the ultrarelativistic outflow is invoked. Aloy 
et al. I|2005p showed by general relativistic calculations 
that the fireball-torus interaction can indeed produce col- 
limated jets with "beaming fractions" /si = 2<5Slj Ct /(47r) = 
1 — cos 6 j C t between ~ 0.5 and ~ 2% of the sky (50j e t is the 
solid angle subtended by one jet and 0j Ct its semi-opening 
angle) . This reduces the true energy outputs of the GRBs 
compared to the observed ones by a factor of 50-200. For 
the recently detected, short, hard bursts GRB 050509b 
(Gehrels et al. 1317051 Bloom et al. 130001 Hjorth et al. l3005jl 
and GRB 050724 (Berger et al. 13005}) this corresponds to 
true 7-energies between some 10 46 ergs and some 10 48 ergs, 
provided the collimation factor of about 100 inferred by 
Berger et al. lUSDSjl for GRB 050724 — which is in agree- 
ment with the predictions of Aloy et al. ( 200]j|) — applies 
to both bursts. Our results suggest that such energies can 
be delivered by vv annihilation without the need to invoke 
additional power from mechanisms based on the presence 
of extremely strong magnetic fields. We therefore refrain 
from discussing here once more possible ways to use such 
magnetic fields for tapping the rotational energy of the 
torus or black hole to drive relativistic motion of small 
amounts of baryonic matter. We feel unable to provide 
any deeper insight into such possibilities on grounds of our 
current simulations than that obtained in the discussions 
of previous works. These exploit the results of hydrody- 
namic models to come up with estimates of equilibration 
field strengths in the swirling gas of the merger remnant, 
or parametrise the possible energy conversion due to the 
Blandford-Znajek mechanism in terms of (undetermined 
or free) system parameters (see, e.g., Rosswog et al.)2003 
Rosswog h Ramirez- Ruiz 1200212003 Lee et al. l2005a| . 
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Besides limited numerical resolution, another draw- 
back of 3D simulations is the very limited time over which 
the evolution of the accretion torus can be followed with 
the computer resources available to us. The models pre- 
sented here certainly add only a small brick to our slowly 
growing picture of the scenario and physics that can make 
up a GRB engine. The ultimate goal certainly must be 
simulations which trace the merger history consistently 
through all phases, from the merging of the compact ob- 
jects, through the formation of the black hole, to the accre- 
tion of the torus matter, until the completion of the phase 
of powerful energy release from this accretion process. 
Such simulations will require general relativistic modelling 
(or a reasonable approximation of it) and should include 
the microphysics which describes the gas properties as well 
as the processes that establish the energy loss of the gas 
and the conversion of some of this energy to relativistic 
and non-relativistic outflows. Further progress in mod- 
elling will not only require adding a spectral treatment of 
the neutrino emission, preferably by solving the 3D trans- 
port problem. It may at some stage also be necessary to 
perform magnetohydrodynamic simulations for investigat- 
ing the combined effects of neutrino transport and mag- 
netic fields. 
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